Ab Initio determination of Cu 3d orbital energies in layered copper oxides 
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It has long been argued that the minimal model to describe the low-energy physics of the high T c superconducting 
cuprates must include copper states of other symmetries besides the canonical 3d x 2_ y 2 one, in particular the 3d z 2 
orbital. Experimental and theoretical estimates of the energy splitting of these states vary widely. With a novel 
ab initio quantum chemical computational scheme we determine these energies for a range of copper-oxides and - 
oxyclorides, determine trends with the apical Cu-ligand distances and hnd excellent agreement with recent Resonant 
"q ' Inelastic X-ray Scattering measurements, available for ha-iCuO^, Sr^CuOnCli, and CaCuC>2. 
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It is generally accepted that the low-energy physics of the layered Cu oxide compounds in their normal state is 
reasonably well described by models which incorporate the "in-plane" Cu 3d x 2_ y 2 and O 2p x /2p y orbitals. However, 
the energy window over which such models provide a qualitatively correct picture is a matter of active research. One 
additional ingredient which is often invoked is the Cu 3d z 2 orbital, perpendicular onto the C11O2 layers, and the apical 
O 2p z functions having c-type overlap with the Cu 3d z 2 . Recent multi-orbital calculations using dynamical mean- field 
theory [l| show indeed that some of the features of the optical, x-ray absorption, and photoemission spectra can be 
better reproduced when the Cu 3d z 2 orbitals are explicitly included in the many-body treatment. At finite doping, 
the inclusion of the Cu 3d z 2 functions makes a difference even for the low-energy states close to the Fermi level [l| . 

The off-diagonal coupling between states of x 2 — y 2 and z 2 symmetry was actually found to substantially affect 
the dispersion of the low-energy bands and the shape of the Fermi surface in earlier semiphenomenological models 
0, 01], density- functional calculations 0), and quantum chemical studies 0. Moreover, Ohta et al. @ and recently 
Sakakibara et al. [7] suggested that a direct relation exists between the magnitude of T c and the size of the d x 2_ y 2- 
d z 2 splitting. The splittings within the Cu 3d shell are also relevant to excitonic models for pairing and high-T c 
, superconductivity [H 0] . Even if the importance of the Cu 3d z 2 state is stressed in this considerable body of work, 
■ the actual experimental and theoretical estimates of the energy of this state vary widely. 

Sharp features at about 0.4 eV in early optical measurements on La2Cu04 and Sr2Cu02Ci2 were initially assigned 
to crystal-field Cu d x 2_ y 2 to d z 2 charge excitations [Toj . A different interpretation in terms of magnetic excitations 
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was proposed by Lorenzana and Sawatzky and latter on confirmed by analysis of the resonant inelastic x-ray 



scattering (RIXS) spectra at the Cu K and L3-edge |12|,|13(. The RIXS experiments also show that in La2Cu04 and 



S^CuC^Ck" the Cu d x 2_ y 2 to d z 2 transitions occur at 1.5-2.0 eV [3], which is substantially larger than the outcome 
of earlier wavefunction-based quantum chemical calculations, 1.0-1.2 eV [HI, or density-functional estimates, 0.9 eV 
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With the aim to settle this point we employ a recently developed ab initio quantum chemical computational 
. £h ! scheme to extract the splittings within the Cu 3d shell in several layered copper oxides. Excellent agreement is 
^ ■ found for La2Cu04, Sr 2 Cu02Cl2, and CaCu02 with recent RIXS measurements Further, the d x 2_ y 2 to d z 2 

excitation energies computed here for La2Cu04, YBa2Cu306, and HgBa2Cu04 are relevant to models which attempt 
to establish a direct relation between the relative energy of the out-of-plane d z 2 level and the critical temperature T c 
0- In particular, the large difference between the critical superconducting temperatures of La2Cu04 and HgBa2Cu04 
was directly attributed to a large difference between the d x 2_ y 2 to d z 2 excitation energies Q. 



RESULTS 



To study bound, excitoniclike states such as the d-d charge excitations in copper oxides, we rely on real-space 
ab initio methods. In the spirit of modern multi-scale electronic-structure approaches, we describe a given region 
around a central Cu site by advanced quantum chemical many-body techniques while the remaining part of the solid 
is modeled at the Hartree-Fock level. The complete-active-space self-consistent-field (CASSCF) method was used to 
generate multireference wavefunctions for further configuration-interaction (CI) calculations [lfj. In the CASSCF 
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scheme, a full CI is carried out within a limited set of "active" orbitals, i.e., all possible occupations are allowed 
for those active orbitals. The active orbital set includes in our study all 3d functions at the central Cu site and 
the 3d x 2_ y 2 functions of the Cu nearest neighbor (NN) ions. Strong correlations among the 3d electrons are thus 
accurately described. The final CI calculations incorporate all single and double excitations from the Cu 3s,3p,3d 
and O 2p orbitals on a given CUO4 plaquette and from the 3d x 2_ y 2 orbitals of the Cu NN's. Such a CI treatment 
is referred to as SDCI. The CASSCF and SDCI investigations were performed with the molpro quantum chemical 



software [17 1 



Both SDCI and RIXS results for the Cu <i-level splittings are listed in Table I. The relative energies of the peaks 
observed between 1 and 3 eV in the Cu L^-edge RIXS spectra [3) are the sum of a crystal-field contribution, 
i.e., an on-site crystal-field splitting E c f, and a magnetic term Ai? mgn . The quantum chemical calculations have 
been performed to extract E c f for a ferromagnetic (FM) arrangement of the Cu d spins. A SDCI treatment for an 
antiferromagnetic (AF) alignment of the Cu spins in the embedded cluster of five Cu sites is computationally not 
feasible (see Methods for details) . AE mgn accounts for AF order in the ground-state configuration of the Heisenberg 
antiferromagnet and is determined as follows. First the value for the NN exchange coupling constant J is computed by 
considering an embedded cluster consisting of two CUO4 plaquettes. For CaCu02, for example, we find J = 0.13eV, in 
good agreement with the theoretical results reported in Ref. [3] and with values from experimental data 11 , 1^, EH ■ 
With this value of J in hand we return to the cluster with five Cu sites and flip the spin of the central Cu ion. 
This corresponds to an energy increase AE — z J/2 — 2J, where z — 4 is the number of NN's and we neglect the 
quantum fluctuations. For the crystal- field excited states, the superexchange with the NN Cu d x 2_ y 2 spins is much 
weaker for a hole excited into the d z 2 orbital and zero by symmetry for a hole into a tig orbital. This contribution 
due to intersite d z 2—d x 2_ y 2 superexchange, AE' = 2J', is not included either in the quantum chemical calculations 
but in a first approximation we can neglect the weak intersite AF interaction J' involving a d z 2 hole. From overlap 
considerations, J' is only a small fraction of the ground-state superexchange J. For a meaningful comparison between 
the SDCI and RIXS data, we subtracted in Table I from the relative RIXS energies reported in Ref. [lj| the term 
AE mgn = AE — AE' w 2 J representing the magnetic stabilization of the ground-state configuration with respect to 
the crystal-field excited states. Since J«0.13 eV, A.E mgn «0.26. 

The agreement between our SDCI excitation energies and the results from RIXS is remarkable. As shown in 
Table I, the differences between the SDCI and RIXS energies are not larger than 0.15 eV. The only exception is the 
splitting between the x 2 —y 2 and xz/yz levels in CaCu02, where the SDCI value is 0.3 eV larger than in the RIXS 
measurements. That an accurate description of neighbors beyond the first ligand coordination shell is crucial is clear 
from the comparison between our and earlier quantum chemical data. In the calculations described in Ref. 15], only 
one CuOe octahedron or one CuOs pyramid was treated at the all-electron level. Farther neighbors were described 
by either atomic model potentials or point charges. Deviations of 0.4 and 0.6 eV (up to 50%) for the z 2 levels in 
La2Cu04 and Sr2Cu02Cl2 [15], for example, are mainly due to such approximations in the modeling of the nearby 
surroundings. The <i-level splittings depend after all on the charge distribution at the NN ligand sites. The latter is 
obviously sensitive to the manner in which other species in the immediate neighborhood are modeled. The quality 
of the results reported here is directly related to the size of the clusters, i.e., five CUO4 plaquettes, all apical ligands 
plus the NN closed-shell metal ions of the central polyehdron. 

Superconductivity has not been observed in Sr2Cu02Cl2 and CaCu02- The d-level splittings for three representative 
cupratc superconductors, i.e., La2Cu04, HgBa2Cu04, and YBa2Cu30g are listed in Table II. The maximum T c 's 
achieved by doping in these three materials are 35, 95, and 50 K, respectively. For the YBa2Cu306 compound, we 
here refer to the maximum T c which can be achieved by Ca doping [l9j. The large difference between the critical 
temperatures in La2Cu04 and HgBa2Cu04 was assigned in Ref. [7| to a large difference between the relative energies 
of the z 2 states in the two materials. The density-functional results for the splittings between the x 2 —y 2 and z 2 levels 
in La2Cu04 and HgBa2Cu04 are 0.91 and 2.19 eV, respectively Q- RIXS data are not available for HgBa2Cu04 and 
independent estimates for the energy separation between the x 2 —y 2 and z 2 states are therefore desirable. While we 
find a rather similar value for HgBa2Cu04, of 2.09 eV, the quantum chemical and RIXS results [3] for La2Cu04 are 
substantially larger, about 1.4 eV. This makes the difference between the d- level splittings in the above mentioned 
compounds less spectacular, i.e., £' H 2 BCO — E L p° is reduced from 1.3 eV in Ref. [7J to 0.7 eV in the present study, 
which suggests that the model constructed and the conclusions drawn in Ref. 0] at least require extra analysis. 
The distance between the Cu and apical ligand sites increases from 2.40 A in La 2 Cu0 4 to 2.78 A in HgBa 2 Cu0 4 



21]. The effect of this growth of the apical Cu-0 bond length on the relative energy of the z 2 hole state can be 



understood by using simple electrostatic arguments: when the negative apical ions are closer to the Cu site, less 
energy is needed to promote the Cu 3d hole into the z 2 orbital pointing toward those apical ligands. For HgBa2Cu04, 
the lowest crystal- field excitation is therefore to the xy level and requires about 1.3 eV, see Table II, while the z 2 and 
xz/yz levels are nearly degenerate and more than 0.5 eV higher in energy. On the other hand, in La2Cu04 the lowest 
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crystal-field excitation is to the z 2 orbital, see Table I. Our results also reproduce the near degeneracy between the z 
and xy levels in La2Cu04, as found in the RIXS experiments. In CaCu02, there are no apical ligands. The splitting 
between the x 2 — y 2 and z 2 levels is therefore the largest for CaCu0 2 , about 2.4 eV, see Table I. 

DISCUSSION 

The parameter that plays the major role in determining the size of the d-level splittings in layered cuprates is 
clearly the apical Cu-ligand distance. There are, however, few other factors which come into play such as the 
number and nature of the apical ligands, the in-plane Cu-0 bond lenghts, buckling of the Cu02 planes, and the 
configuration of the farther surroundings. Trends concerning the relative energy of the z 2 hole state in different 
cuprates are illustrated in Fig. 1, which includes data for systems having one apical O site (YBa2Cu30g), two apical 
O's (La2Cu04, HgBa2CuC>4), two apical CI ions (Ca 2 Cu02Cl2, S^CuC^C^) or no apical ligand (CaCuO?). The 
apical Cu-0 distances in La2Cu04 and YBa2Cu306, for example, are nearly the same, 2.40 vs. 2.45 A [lilliol 22|. 



In YBa2Cu306, however, there is a single apical O. For this reason the z 2 hole state is somewhat destabilized in 
YBa2Cu306 and lies above the xy hole configuration, see Table II. Yet since the Cu ion is shifted towards the apical 
ion, out of the basal O plane, the x 2 —y 2 hole state is also destabilized such that the splitting between the x 2 —y 2 and 
z 2 levels is finally close to the value found in La2Cu04. Further, the apical Cu-ligand distances are slightly larger in 
Sr2Cu02Cl2 as compared to HgBa2Cu04, 2.86 vs. 2.78 A, respectively. The apical ions also have a smaller effective 
charge in Sr2Cu02Cl2, which should lead to a larger relative energy of the z 2 hole state in Sr2Cu02Ci2 as compared to 
HgBa2Cu04. The fact that the relative energy of the z 2 hole state is actually larger in HgBa2Cu04, see Fig. 1, must 
be related to the smaller in-plane Cu-0 distances in HgBa2Cu04, 1.94 in HgBa2Cu04 vs. 1.99 A in Sr2Cu02Cl2, 
which stabilizes the ground-state x 2 —y 2 hole configuration in the former compound, and to farther structural details. 



From Ca2Cu02Cl2 to Sr2Cu02Cl2, the Cu-Cl separation increases from 2.75 to 2.86 A 23j, |24| and the energy of the 
z 2 level from 1.37 to 1.75 eV. 

In contrast to the z 2 orbitals, the relative energies of the xy levels display much smaller variations, in an interval 
of 1.2—1.5 eV, see Tables I and II. Substantially smaller are also the variations computed for the xz/yz levels, in an 
energy window between 1.6 and 2.0 eV. 

To summarize, we employ state of the art quantum chemical methods to investigate the Cu 3d electronic structure 
of layered Cu oxides. Multiconfiguration and multireference configuration-interaction calculations are carried out 
on finite clusters including five CU.O4 plaquettes plus additional apical ligand and closed-shell metal ion NN's. The 
localized Wannier functions attached to these atomic sites are obtained from prior Hartree-Fock computations for the 
periodic system. Excellent agreement is found between our theoretical results and recent Cu L3-edge RIXS data for 
La2Cu04, Sr2Cu02Cl2, and CaCu02. RIXS is a novel experimental tool to investigate both magnetic and charge 
excitations with high resolution and accuracy. Our computational scheme and present results indicate a promising 
route for the modeling and reliable interpretation of RIXS spectra in correlated 3d-metal compounds. A next step 
along this path is the computation of transition probabilities and intensities at the ab initio level, which requires the 
explicit calculation of the intermediate Cu 3p core hole wavefunctions. 

Further, the excitation energies computed here for La2Cu04, YBa2Cu306, and HgBa2Cu04 are relevant to models 
which attempt to establish a direct relation between the critical temperature T c and the strength of the (x 2 —y 2 ) — z 2 
coupling. For La2Cu04, in particular, the density-functional estimate used as input parameter in such models Q 
is about 0.5 eV smaller than our result. Consequently, the difference we find between the d x 2_ y 2-d z 2 splittings in 
La2Cu04 and HgBa 2 Cu04 is less spectacular as compared to the value reported by Sakakibara et al. [7], suggesting 
a reevaluation of the analysis in Ref . [7[ . 

METHODS 

The first step in our study is a restricted Hartree-Fock (RHF) calculation for the ground-state configuration of the 
periodic system. The RHF calculations are performed with the crystal package [25j. We employed experimental 
lattice parameters 14, 23-24| and Gaussian- type atomic basis sets, i.e., triple-zeta basis sets from the crystal 



library for Cu, O, and CI plus basis sets of either double-zeta or triple-zeta quality for the other species. Post Hartree- 
Fock many-body calculations are subsequently carried out on finite clusters, which are sufficient because of the local 
character of the correlation hole. They consist of five CUO4 plaquettes, i.e., a "central" CUO4 unit plus the four 
NN plaquettes. When present, the apical ligands, oxygen or chlorine, are incorporated as well in the finite cluster C. 
Additionally, the finite cluster C includes in each case the NN closed-shell metal ions around the "central" Cu site. In 
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La2Cu04, for example, there are ten La 3+ NN's. In YBa 2 Cu 3 6 , there are one Cu 1+ 3d 10 , four Y 3+ , and four Ba 2+ 
NN's. 

The orbital basis entering the post Hartree-Fock correlation treatment is a set of projected RHF Wannier functions: 
localized Wannier orbitals (WO's) are first obtained with the Wannier-Boys localization module (2(| of the CRYSTAL 
package and subsequently projected onto the set of Gaussian basis functions associated with the atomic sites of C [27 1. 
Moreover, the RHF data is used to generate an effective embedding potential for the five-plaquette fragment C. This 
potential is obtained from the Fock operator in the RHF calculation j2?| and models the surroundings of the finite 
cluster, i.e., the remaining of the crystalline lattice. 

The central CUO4 plaquette and the four NN Cu sites form the active region of the cluster, which we denote as 
Ca- The other ions in C, i.e., each ligand coordination cage around the four Cu NN's and the NN closed-shell metal 
ions, form a buffer region Cb whose role is to ensure an accurate description of the tails of the WO's centered in the 



active part Ca [27J. For our choice of Cb, the norms of the projected WO's centered within the active region Ca are 



not lower than 99.5% of the original crystal WO's. While the occupied WO's in the buffer zone are kept frozen, all 
valence orbitals centered at O and Cu sites in Ca (and their tails in Cb) are further reoptimized in multiconfiguration 
CASSCF calculations. In the latter, the ground-state wavefunction and the lowest four crystal-field excited states at 
the central Cu site are computed simultaneously in a state-averaged multiroot calculation [l6| . The d-level splittings 
at the central Cu site are finally obtained at the CASSCF+SDCI level of theory as the relative energies of the crystal- 
field excited states. The virtual orbital space in the multireference SDCI calculation cannot be presently restricted 
just to the Ca region. It thus includes virtual orbitals in both Ca and Cb, which leads to very large SDCI expansions, 
~ 10 9 Slater determinants for a FM configuration. For this reason, we restrict the CASSCF+SDCI calculations to 
FM allignment of the Cu d spins. 

The effective embedding potential is added to the one-electron Hamiltonian with the help of the CRYSTAL-molpro 
interface program j28|. Although the WO's at the atomic sites of C are derived for each of the compounds discussed 
here by periodic RHF calculations for the Cu 3d 9 electron configuration, the embedding potentials are obtained by 
replacing the Cu 2+ 3c? 9 ions by closed-shell Zn 2+ 3g? 10 species. This is a good approximation for the farther 3d-metal 
sites, as the comparison between our results and RIXS data shows. An extension of our embedding scheme toward 
the construction of open-shell embcddings is planned for the near future. 
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FIGURE LEGENDS 

Figure 1. Relative energy of the Cu z 2 level as function of the distance h between the Cu and apical ligands in 
different cuprates. There are two apical O sites in La2Cu04 and HgBa2Cu04, two apical CI sites in Ca2Cu02Cl2 and 
Sr 2 Cu02Cl2, and one apical O ion in YBa2Cu 3 06- In CaCuC-2 the apical ligands are absent. 
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TABLE I: CASSCF+SDCI versus RIXS results for the Cu d-level splittings in La 2 Cu0 4 , Sr 2 Cu0 2 Cl2, and CaCu0 2 (eV). 
The ground-state Cu fynd^d^z-y" 1 configuration is taken as reference. A 2 J term was here substracted from each of the RIXS 
values reported in Ref. [lifl. see text. 
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